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Rheology of a granular gas under a plane shear 
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The rheology of a two-dimensional granular gas under a plane shear is investigated. From the 
comparison among the discrete element method, the simulation of a set of hydrodynamic equation, 
and the analytic solution of the steady equation of the hydrodynamic equations, it is confirmed that 
the fluid equations derived from the kinetic theory give us accurate results even in relatively high 



I. INTRODUCTION 

Granular materials consist of macroscopic dissipative particles. In some cases the granular material behaves as an 
i-Q ' unusual fluid. Although to understand the rheology of the granular fluid is practically important, our understanding 
Ti \ on the rheology is not enough. There are several reasons to have poor understanding on the rheology of granular flows: 
(i) The separation of the length scale between particles and the fluid motion is not enough, (ii) there are some cases 
that the fluid region coexists with the solid-like region, and (iii) most of experiments are strongly affected by boundary 
conditions and the external field. Nevertheless, it is believed that rapid granular flows for relatively dilute granular 
gases can be described by a set of hydrodynamic equations at Navier-Stokes order whose transport coefficients can 



nyc 
^ . be calculated by the kinetic theory. Q 



To maintain a granular gas we need to add an external field. The simplest steady situation of the granular fluid is 

achieved by the balance between an external shear and inelastic collisions between particles. This system is appropriate 

C ■ to investigate what the constitutive equation for the granular fluid is. The kinetic theory may suggest that the stress- 

i-rt strain relation can be described by that at Navier-Stokes order, though the transport coefficients can be functions of 

^ the position. 10,11] 

O About 50 years ago, Bagnold[3 suggested that the granular fluid is characterized by r ex 7^ where r and 7 are 

, O, . the shear stress and the shear rate (the strain) , respectively. This stress-stain relation is known as Bagnold's scaling 

and is different from the conventional Newtonian relation which is r ex 7. Recently, Pouliquen|y| and Silbert et 

P\J , al.\j^ have reconfirmed the quantitative relevancy of Bagnold's scaling in granular flows on inclined slopes. Mitarai 

K*' ' and Nakanishijsl have demonstrated that the kinetic theory can be compatible with Bagnold's scaling, when they 

^^ ] assume that the the temperature is a slaving variable of the velocity and the density. Santos et a/.j^] also indicate 

that Bagnold's scaling is valid for steady dilute granular gases without the influence of the gravity in the uniform 

shear flow (USF), though the transport coefficients such as the viscosity and the heat conductivity are different from 

^~j . those in homogeneous cooling stateQ. However, we still do not know whether Bagnold's scaling is relevant in other 

\^ ' situations. 

f^ , We should recall that it is difficult to keep granular gases in experimentally relevant situations because of the 

^^ • existence of gravity. Recently, to remove the effects of gravity, the rheology of dense granular flows under the 

C^ [ plane shear with a constant pressure has been studied and a new scaling has been reported [ifl [13 ■ These studies are 

H I important but particles are not in a gas state, i.e. each particle is in contact with many other particles simultaneously. 

I The analysis of such process is challenging but we do not have any good tool to analyze it at present. Here, we focus on 

the granular shear flows without multi-body contacts in a constant volume container to discuss quantitative relevancy 

of the kinetic theory. 

^ , The purpose of this paper is to know the relevancy of the kinetic theory for a granular gas with moderate density 

under the plane shear in a constant volume container to characterize the rheology of granular fluids. For this purpose, 

we investigate whether (i) the kinetic theory is relevant for the relatively dense granular gas, and (ii) the tangential 

contact and the rotation of particles are irrelevant except for the boundary layers. It should be noted that Bagnold's 

5h scaling is no more relevant for flows with non-uniform shear rate in the steady granular flow. In fact, to satisfy r ex 7^, 

the shear rate 7 should be uniform, because the shear stress and the pressure are spatially uniform in the steady state. 

To investigate the above problem, we use the discrete element method (DEM) for particles' simulation in a plane shear 

problem of dense granular gases in which particles are conflned in a constant volume box (Section II). We adopt the 
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constitutive equation of hydrodynamics derived by Jenkins and Richman 12^ for non-rotational particles in Section III. 
In Section IV we solve a set of hydrodynamic equations with the above constitutive equation numerically and compare 
the result with the result of DEM. We also check the renormalization theory of the restitution coefficient developed 
by Yoon and Jenkins 13]. In Section V we obtain the analytic solution of the steady hydrodynamic equations to 
verify the quantitative relevancy of the constitutive equation. In section VI we discuss our result and the relevancy of 
Bagnold's scaling. In Section VII, we conclude our results. In Appendix, we briefly explain the method to determine 
the tangential restitution constant as a function of the incident angle. 

II. DEM SIMULATION 
A. DEM model 

The discrete element method (DEM) is one of the standard methods to simulate granular motions. [l4J DEM is 
applicable to most of situations of granular dynamics even when particles are almost motionless and contact with 
many other particles. We adopt DEM to simulate a granular fluid to check (i) the validity of kinetic theory based on 
the reliable model, and (ii) the effects of rotation of particles for the granular fluid. 

In this paper, we focus on a two-dimensional motion of granular particles under a plane shear. We adopt the linear 
spring model for the repulsion with the normal stiffness fc„ and the tangential stiffness kt, and the normal and the 
tangential viscous coefficients ?7„ and ?]t, respectively. 

Let us consider a colliding pair of two disks i and j of the diameter a and the mass m at the position x^ with the 
velocity Ck = Xk and the angular velocity tjJk for k = i or j. If the particles are in contact, the overlap distance 

A,^ EE 2a ~ \x, - Xj\ (1) 

must be positive. The relative velocity at the contact point is 



-n.y X (uji +u;j) (2) 



with the normal unit vector n^ = {xi — Xj)/\xi — Xj\. Introducing the normal velocity cjf = n^ • Cij, the tangential 
velocity c^ = tij ■ Cij, the tangential displacement Wi = /. dsci {s) with the tangential unit vector i^ satisfying 
tij ■ riij — 0, the normal and the tangential forces F^^ and Fj:-' are respectively given by 

F^^ = r7ifc„Ay - m?7„u,Y for Ay > 0, (3) 

F^' = 7mniht,fi\F^^)stgn{hl^), (4) 

where h]'' = —mktwl^ — mrjic^ with Coulomb friction constant /i, min{a, b) is the function to select the smaller one 
between a and 5, and sign{x) — 1 for a; > and sign{x) = — 1 for x < 0. The total repulsive force at the contact can 
be represented as Fij = F^^riij + F^Hij. 

Thus, the equation of motion of particle i is described by 



mci 



E^^^' (5) 



j#i 



^"^^ " |E"y ^*y^t'^ (6) 

where I = ma"^ /8 is the moment of inertia. We integrate (jS)) and @ in terms of the second order Adams-Bashforth 
with the time interval 6t = 4.0 x 10^'^{2a/U). 

Through the paper we adopt the following parameters as /c„ = 3.0 x 10^{U/a)^, kt = fc„/4, 7y„ = 3.0{U/a), i]t = '7n/2 
and II = 0.20, where U is the shear speed at the boundary. These parameters lead to the normal restitution constant 
e = 0.85 and the tangential restitution /3 ~ —1 -I- 1.12442 cot7 for 7 < 7c and /3 = /So — 0.769235 for 7 > 7c where 
7 is the incident angle of two colliding disks and the critical angle 7c is given by cot 7c ~ 1.56734 (see ea. (|A4|l in 
Appendix). As shown in Appendix, the tangential restitution constant (3 can be approximated bv(l5| 

g^i-l + Kl + e)cotj{l + ^) (7>7c) ,^. 

\po (7 < 7c). 

We also note that the realistic value of Coulomb friction constant in both disks and spheres is /x < 0.2 [1^, u^ u^ 
Thus, the renormalization theory of the restitution constant may be applicable to many realistic situations. 



B. Set-up 
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FIG. 1: The configuration of our simulation. 



Simulations of granular particles under the plane shear have been performed by many researchers, but many of 
them[l^ I23, Um assume Lees-Edwards boundary condition[23| which may not be adequate to consider the effects 
of physical boundary. On the other hand, Babic[23|, Popken and ClearyJ2J| have simulated sheared granular flows 
confined in frictional flat boundaries, but their simulations are restricted to the cases for small systems and almost 
elastic particles. Kim|23 has indicated that the density of particles near the boundary is higher than that in the 
bulk region for his simulation of a small system with the flat frictional boundary, while particles are accumulated in 
the center region for a larger system. Louge 26] has simulated a three dimensional shear flow on the flat frictional 
boundary to examine the boundary condition proposed by Jenkins |27l|. but Louge is mainly interested in the behavior 
of flux, the stress ratio as the functions of volume fraction and the restitution constant. Recent papers by Xu et 
al. for an ex peri ment |28j| and a simulation J23 examine the validity of three dimensional kinetic theory by Jenkins 
and Richman[30j under asymmetric shears in the presence of a streamwise body force, where they obtain reasonable 
agreements between the theory and the observations in both the experiment and the simulation. 

As long as our knowledge, we do not know papers to discuss the validity of kinetic theory in a transient dynamics and 
a symmetric shear without a streamwise body force with an enough large system. Thus, we adopt the following setup 
of our DEM simulation shown in Fig^ The system is confined in a two-dimensional container. Without including 
the effects of the air and the gravity we add a symmetric shear with the shear speed U as shown in Fig^ The 
parameters are fixed as the number of particles N = 5000, the linear dimension of the system in y direction A = ISOcr 
and the mean area fraction v — 0.I2I. The boundary condition of x— direction is periodic. We non-dimcnsionalizc all 
quantities by the diameter a for the length scale, m for the mass, and the inverse of the shear rate 2a /U for the time 
scale. 

We introduce some fixed particles on the wall to reproduce the bumpy boundary. The reason why we adopt the 
bumpy boundary is to avoid the large amount of slips on the wall under a physical situation. In our simulation we 
start from an initial condition without the shear. Then, the wall at y = A/2 obeys the equation of motion in x 
direction 



M, 



dcw 

'~dr 



-mj^{cw - Ue^/2) + F^ 



(8) 



where cw and e^ are respectively the actual wall velocity and the unit vector along a;— direction. M^, 7^ and F^x 
are the mass of wall M^, = 5.0 x IO^tti, the relaxation rate 7^, = 10C//(2ct), and the force acting on the wall by the 
collision between mobile particles and the wall, respectively. 



Simulation 





(a) (b) 

FIG. 2: The time evolution of tiie particles' configurations from (a) to (b) for v = 0.121. 

The initial condition is prepared as that the configuration of particles is at random and the velocity distribution 
function obeys Maxwellian. Figure |21 is the time evolution of particles' configuration for v = 0.121. Two shallow 
clusters appear near the wall, and move to the center region of the container. Then, the two clusters merge to form 
a big cluster. A similar behavior can be observed in the simulation under the Lees-Edwards boundary condition. |2l| 
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FIG. 3: The time evolution of the kinetic energy, where the units of time and the energy are 2a jU and mU'^ /2, respectively. 



by 



The time evolution of the total kinetic energy E m. a, typical situation is shown in Fig|31 where the energy is defined 



E{t)^\Y,{m^,+lLol). 



(9) 



In this figure, the initial energy is larger than the steady value, but we have checked that the qualitatively similar 
results can be obtained even when we start from the smaller energy about E{G) = 1200 in the dimensionless unit. It 
is characteristics that the total kinetic energy is relaxed to be almost a constant value quickly, but there is the slow 
evolution of hydrodynamic fields. 

The hydrodynamic variables are the local area fraction v(r,t) = 7ro'^n(r, t)/4 with the number density n{r,t), the 



velocity field v{r, t) and the granular temperature T{r, t) which is defined by 

T{r,t)^^Jdc{c-vff{r,c,t). (10) 

In our simulation we divide the system into square cells with the linear dimension a. Then, we can check to what cell 
each particle belongs. Thus, the measured area fraction in our simulation is given by i^(r,i) = X^iec "''-'^^/^ where 
the summation is taken over the center of particle i existing in the cell C at r with the area A — a^. Similarly, the 
velocity field v{r, t) = Yln^c ^i ^^ ^^^ local average of the velocity of particles. The temparture field is also calculated 
by 

1 



Tir.t) 



2n ^ 



(c. 



(11) 
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FIG. 4: The time evolution of dimensionless Uy with (a) at f = 60 and (b) at f = 380, where the units of the time and Uy are 

{2a/U) and C//2, respectively. 

The time evolution of hydrodynamic variables can be summarized as follows. Corresponding to Fig. 2, the local area 
fraction becomes large near the boundary at the initial stage, and the shallow clusters move to the central region of 
the container with growing the peak density. Finally, the two clusters merge to form a compact cluster at the center 
y = 0. 

The interesting behavior can be observed in the velocity field and the granular temperature. The x-component of 
the velocity field is almost zero between two clusters during the time evolution, and such plug region is narrower as 
the distance of two clusters is closer. Even in the steady state, the velocity gradient in the central region is smaller 
than that near the boundary. The behavior of y-component of the dimensionless velocity field Uy is also characteristic, 
because this quantity shows the tendency to move to the central region and is relaxed to be zero in the steady state 
(Fig^J. Similarly, the temperature field in the central region is smaller than that near the boundary. However, the 
minimum value of temperature decreases with time and reaches almost zero in the central region in the steady stage. 
The results of our DEM simulation except for Uy — 2vy/U will be shown in Sections IV and V through the comparison 
of DEM with hydrodynamic simulations or the theoretical results in the steady state. 

The result of our simulation may give us suspicious impression of the validity of kinetic theory for this system, 
because (i) the density in the cluster is near the closest packing 7r/2-\/3 ~ 0.907 in the steady state, (u) the particles 
are almost motionless in the clusters or between clusters. Our result is contrast to the result by Xu et aL|28l l29l | 
where there is no definite motionless clusters because of the existence of a streamwise body force. In their case, the 
mean density is much higher than that in our case, and the temperature is a nearly constant. Nevertheless, as will 
be shown, our system can be described by the hydrodynamic equation derived from the kinetic theory. 



III. HYDRODYNAMIC EQUATIONS 



The purpose of this paper is to verify the validity of the granular hydrodynamic equation derived from the kinetic 
theory. Although there are two standard methods, Chapman-Enskog method 0,|jj|31| and Grad expansions 32], most 



of established results are limited to low density cases in the derivation of hydrodynamic equations. However, as shown 
in the previous section, we have to adopt the kinetic theory for dense granular gases. The Chapman- Enskog method 
by Grazo and Duftyj^ and Lutsko[23 predicts transport coefficients in dense granular gases. On the other hand, 
Jenkins and Richman derive hydrodynamic equations based on Grad expansion and give transport coefficients 12, SO"]. 
Although the treatment by Jenkins and Richman JJj does not take into account the contraction of the phase space 
volume in each collision, the theory is suitable for our purpose because it gives us explicit expressions of the transport 
coefficients in the two-dimensional dense granular gases. 

Jenkins and Richmanil^], and Lun 35] derive sets of hydrodynamic equations which include the angular velocity, 
the spin temperature as well as the density, the translational velocity and the granular temperature. The equations 
include the couple stress and the collisional loss of spin energy. These sets of hydrodynamic equations are categorized 
as the micropolar fluid mechanics which was originally proposed by Cosserat and Cosserat for the description of 
the elastic materials |36j. Application of the concept of micropolar fluid mechanics to atomic gases are developed by 
Dahler and his coworkers |3^. The micropolar fluid mechanics is applied to granular flows by KanatanipJSJ. Lun|35|. 
Babic 39|, HayakawaW^ and Mitarai et aZ.|4^. The importance of the excitation of the spin on the boundary is 
indicated by Jenkins 27], but the effects of the spin can be decoupled with the translational velocity in the bulk 
region. In particular, Babic^39j indicated that the coupled stress induced by the collisions between circular particles is 
canceled. Therefore, we believe that the effect of spins can be absorbed in the boundary condition and the renormalized 
restitution constant. 

Recently, Jenkins and Zhang |42J have proposed a scheme of the renormalized restitution constant as the result of 
the absorption of the rotation of particles. Yoon and Jenkins[l3 have extend the scheme to two-dimensional flows as 

e~e-iJ. + 2^?[l + e). (12) 

The validity of three dimensional theory [43 has been tested by Xu et al.^^ and Jenkins and Zhang p^. The latter 
is consistent with Lun and Bent|43l| in part. However, the quantitative validity of Yoon and Jenkins[l^ has not 
been confirmed yet. In addition^ some recent papers suggest that the spin effects are relevant in granular flows. 
For instance, Goldhirsch et aZ.jij] have indicated that the equipartition between spin energy and the translational 
energy is violated in a recent paper, and Gefen and Alam 45] discuss the linear stability of sheared micropolar fluid. 
Therefore, we need to judge whether the concept of micropolar fluid is necessary for the description of granular fluids. 

In this paper, we adopt the renormalization procedure of the restitution constant proposed by Yoon and Jenkins|l3J 
to verify the validity of their scheme. Thus, the restitution constant e appears in hydrodynamic equations is different 
from e of DEM, where the relation between two restitution constants is given in ea. l|12|) . Namely, e = 0.85 in DEM 
corresponds to e = 0.798 for ^ = 0.20 in hydrodynamic equations. 

The advantage to adopt the renormalization is that hydrodynamic equations can be simplified as 

Dtp = -pV-v, (13) 

pDtV = -V-P, (14) 

pDtT = -P:iVv)-V-q-x, (15) 

where p = nm is the mass density and Dt = dt + v -V . Here {i, j) component Pij of the pressure tensor P is expressed 
as the function of the bulk viscosity ^ and the shear viscosity 77 

P^J ^ [p ~ ^{V ■ v)]S,, ~ ijD,, (16) 

at the Navier-Stokes order, where 6ij = I ioi i — j and for otherwise, Dij — {WiVj + WjVi)/2. Here q represents the 
heat flux which can be expanded as 

q = -kVT - XV p, (17) 

where n is the heat conductivity and the transport coefficient A disappear at e = 1 . The collisional loss rate of the 
energy x can be represented by 

X = -^^^p^ff(^)Ti/2[8T - 3V^aT'/\V ■ v)], (18) 

where pp — Am/ [1:0'^) is the mass density of a particle. 

Let us non-dimensionalize the time, the position, the velocity and the temperature as 

i=^i*, x^ax\ ^^ju^ ^^\^ (^^) 



TABLE I: The dimensionless transport coefficient by Jenkins and Richman 



pM = ^''[^ + i^ + (^>9{'') 



CM = -^{l + eyg{u) 



/tt. 1 , ,_i (l + e)(3e+l) ,(l + e)(3c-l) 1,,^ x 2 , m 

''"> = V 2 It^"*"' + 4(7- 3e) ■' + ' 8(7 -3^) + ;;«' + =>" "("M 
»M ^ ^iiTT^i(kii:j»<'')-'-lS^-'^^fe;i^-i^>('^ 



Thus, the non-dimensional pressure tensor, the heat flux and the collisional loss rate of energy are respectively given 
by 



4 '^' ^ 8 ^ ' ^ 8a 



P^^-^^Ptv 9 = ^9*, X=^X*. (20) 



Here the dimensionless quantities are written as 

P*^ = W)9^i{u)e^'W*-u)]6,,-r^{i.)e^'^b:^, (21) 

q* ^ -K{v)0^'^V*e - \{v)0'^/^V*v, (22) 

X* = ^^2<?(^)0i/2[40-3,&/2(V*.«)]. (23) 



4V2^ V 2 

The explicit expressions oi p{v), S,{v), ri(i/), kJi/) and A(i^) obtained by Jenkins and Richmanll2l| are summarized in 
table m with the radial distribution function[4y] 

.(^)=.c(.) + — ^4^^p^^, (24) 

l + exp[-(t/- t/o)/"ioJ 



where gc(!^) = (l - 7iy / 16) / {1 -1^)^ and gf{i^) = [(1 + 6)1^(^^^-1)]"^ with i/^ = 0.82, i/q = 0.7006 and mo = 0.0111. 
The choice of 5(1^) is not be unique. For example, we expect that a similar result can be obtained by using the radial 
distribution function in ref. 47]. Thus, the dimensionless hydrodynamic equations are reduced to 



Dtv = -vV -u, (25) 

vDtU = -V-P, (26) 

1 
2' 



-yDtO = -P,,{V,Uj)-V-q-x- (27) 



From here, the mark * to represent dimensionless quantities is eliminated. 

IV. SIMULATION OF HYDRODYNAMIC EQUATIONS 
A. The outline of our simulation 

To verify the accuracy of a set of hydrodynamic equations (|25|l - (|27|l derived from the kinetic theory by Jenkins and 



: ity t 

m 



Richman 1 12|. we simulate hydrodynamic equations. Since the separation of particles' scale and the hydrodynamic scale 



8 

is not enough, each grid in two-dimensional space cannot contain enough number of particles to define hydrodynamic 
variables. Therefore, we focus on the field equations which have translational symmetry in a;-direction. Thus, all 
quantities only depend on y and t. However, we should note that we keep x-component of velocity. The second 
purpose of the simulation of hydrodynamic equations is to obtain a reduced set of equations to recover the qualitative 
accurate results to describe the metastable dynamics after the total energy is relaxed to a constant value. 

The method of the discretization of continuous variables is based on the standard procedure. We adopt the classical 
Runge-Kutta scheme for the time derivative with 5t = 0.01 and the second order accuracy of the spatial derivative of 
a hydrodynamic variable \1/ as 



9* _ *j+i - *j_i 92^ _ fj+i - 2*j + * 



j-i 



dy 2h ' dy^ h^ 



(28) 



where h is the grid displacement with h/A = 1/180 and y = jh with j ^ 0, ±1, ±2, • • • . It should be noted that we do 
not have to solve Poisson equation for the pressure because the fluid is compressible and the pressure is completely 
determined by the equation of state. 

B. The boundary condition 

We adopt the boundary condition proposed by Johnson and Jackson 48] . We define the slip velocity on the boundary 
as Usi = u — Uw, where u^ = ie^ at y = ±A/(2cr). Let t and n be the tangential unit vector and the normal unit 
vector to the wall, respectively. Thus, the conservation of the linear momentum on the wall is given by 

-n-P-t=^(j)fliiy,e)\usi\, (29) 

where 7r/4 is originated from m — irppa^ /4. Here, </) is the roughness parameter and i^{v, 0) is the collisional frequency 
between the wall and the particles. The expression fl is assumed to be 

r!(^,0) = ^.5(j.)0i/2, (30) 

where the prefactor is absorbed in the roughness parameter. On the other hand, the energy balance on the wall can 
be expressed by 

n-q=-u,i-P-n-T{iy,6) (31) 

where T^v, 9) is the energy loss rate in terms of the inelastic collisions between particles and the wall, which may be 
represented as 

T{v, 9) = ^<Pn{,y, 6)9 = ^^vg{v)9'"\ (32) 

where $ is the hardness parameter of the wall. In our simulation we adopt (f) — 0.20 and $ = 0.24 as fitting parameters. 

The reason why we adopt the boundary condition by Johnson and Jackson|43 is that their condition is simple. 
The more precise treatment for the boundary condition can be seen in ref.|49|. When we adopt Jenkins' boundary 
condition, the number of fitting parameters may be reduced. 

When we represent these boundary conditions as 

Fw(*,5y*) = (33) 

and the formal solution of this discrete equation can be formally solved as 

*7V = i^b2(*A'-l, *iV-2), (34) 

where N = A/2ct is the grid number on the boundary with symbolic functions Fj^i and Fi,2- From the consideration 
of the symmetry in y-direction, we have v{y) — v{^y), 9{y) = 9{—y), and u{y) — ~u{—y). Thus, it is enough to 
discuss the boundary condition at y = A/2. From 129|l we may obtain the x— component of the velocity field 

Ux;N~2 + ^4)hVN-ig{vN-l)/-n{vN-l) .„_. 

"^;JV = , , , ,, -, ■—; ^ , (35) 

1 + A^hvN-igyvN-ijriivN-i) 

where the area fraction vjq on the boundary is assumed to be 

UN = 2j/Ar_i + l/Ar_2 (36) 
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FIG. 5: The initial conditions of hydrodynamic equations (solid lines) and the corresponding data of DEM (open circles) at 
t = 20. The solid lines for v and Q are the polynomials of even powers of y until y , while the lines for Ux 
polynomials of the odd powers of y until y^. 



to suppress the gradient of the density field. Similarly, ea. (|31|) becomes 



^N 



IN 



-i\{vn-i){vn-2 - Vm) 



2h 



[0(1 - u^-n)'^ 



In 



-i\t^N-ig{vN-i) 



k{vn-i) 
It is obvious that y component of the velocity field satisfies 

Uy;N = 0. 



k{vn-i) 



(37) 



(38) 



C. The result of numerical simulation for the complete set of hydrodynamic equations 

For the initial condition to simulate hydrodynamic equations we fit the data of DEM at t = 20 in the dimensionless 
unit. Each fitting curve is approximated by a polynomial of y (FigEJ. The reason why we adopt the initial condition 
at i = 20 instead of i = 0, we are interested in the slow evolution of hydrodynamic variables after the total kinetic 
energy is relaxed to be a constant. 

As shown in Figs. IHl and [Tithe results of the simulation of hydrodynamic equations well agree with those of DEM. 
Once we rescale the time, the evolution of hydrodynamic variables in the simulation of hydrodynamic equations 
is almost equivalent to that of DEM. This agreement between hydrodynamic equations and DEM means that the 
renormalization procedure by Yoon and Jenkins J13l| gives us the accurate results. Amongst hydrodynamic variables 
the y— component of the velocity field in the simulation of hydrodynamic equations has much larger than that of DEM 
though the profile itself is similar with each other, but the other variables in hydrodynamic simulation are almost the 
same as those in DEM (Fig[7|). 
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FIG. 6: The comparison of the data for the area fraction for v = 0.121 obtained by DEM (open circles) at f = 20 (the label 1), 
60 (the label 3) and 380 (the label 5), and the result of hydrodynamic equations 125ll - 127|l (solid lines). 
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FIG. 7: The time evolution of the granular temperature (a), the velocity fields u^ (b) and Uy (c) in hydrodynamic simulations 
shown in solid lines, where the data in (a) and (b) are obtained by DEM. The numbers 1,3,5 in these figures correspond to 
results at i = 20, 60 and 380, respectively. The DEM data with the solid squares and open circles correspond to the result at 
t = 20 and 60, respectively. Note that comparison of the theory and DEM in the steady values of Q and u^ will be shown in 
the next section, while we do not include DEM data for Uy because of the disagreement in the scale (see Fig. 4). 

It should be noted that the transient dynamics of a granular shear flow has been discussed by Babic[2^ but his 
system is relaxed to be an USF because of the small system size and inelasticity. On the other hand, ours will not 
reach USF, and the time evolution of hydrodynamic variables contains a pattern formation. 



D. Simulation of a simplified set of equations 



To understand the qualitative behavior of phase separations, we need to reduce the degree of freedom of hydrody- 
namic equations. It is reasonable to deduce the terms proportional to the bulk viscosity is not important. In addition, 
the advection term u ■ Vm in hydrodynamic equations may not play important roles because we are interested in the 
slow dynamics in the domain growth. The coupling between the spatial gradient and the terms proportional to 1 — e^ 
because the kinetic theory can be applied to cases for small inelasticity. 



Thus, we may reduce the set of hydrodynamic equations to be 



dtiy 
dtUx 

dtUy 
dtO 



= -dy{vUy), 



= dy[ 



^^\-l-d.,u. 



dy[f]{iy)e'/^dyUy-p{iy)9], 

-UydyO - V-^'q{v)e^/^{dyU^f + 2^/- ^ 7? ( l^) ^ ^ Z^ (^y Uy ) ^ - 2i/ " ^ Sj, ( K ( i^) 6* ^ ^^ ^^ 

-V2^{l-e^)yg{v)e'^'^. 
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(39) 
(40) 
(41) 

(42) 



For further simpUfication, we may assume that the temperature is a fast variable to slave other slow variables. 
Thus, we may omit the time derivative in ea. (|42|l . but such simplification does not lead to a simplified treatment 
to solve hydrodynamic equations. Although Uy becomes zero in the steady state, it plays an important role for the 
time evolution of density. Thus, we believe that the set of equations (|39|I - H42|I is the simplest set of hydrodynamic 
equations to describe phase separations. 

Figure IHl shows the growth of hydrodynamic variables based on H39|l - (|42|l . Although the quantitative behavior is 
a little deviated from the result of DEM or the full set of hydrodynamic equations (|25|I - H27|) in particular for u^, 
qualitative behavior of this simplified model is similar to those in more accurate treatments. In the steady state both 
hydrodynamic models reduce to equivalent results. 
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FIG. 8: The time evolution of the the area fraction (a) granular temperature (b), the velocity fields Ux (c) and Uy (d) in the 
simulation of a simplified model 1421 of hydrodynamic equations. The numbers 1,3,5 in these figures correspond to the data at 
t = 20, 60 and 380, respectively. 
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V. THEORETICAL DESCRIPTION OF THE STEADY STATE 

In the steady state, the hydrodynamic variables depend only on y. Thus, the variables are 

v = v{y), Ux = Ux{y), Uy = 0, e^Oijj). (43) 

It is obvious that any hydrodynamic variable ^ satisfies Dt^ — dt'i = 0. Thus, the equation of mass conservation is 
automatically satisfied. The remain equations of motion become 



d_ 
dy' 
d 
dy' 

d 



= —P^y, (44) 

= ^Pyy, (45) 



Thus, the normal stress and the shear stress are uniform 

p = Pyy — const. T = Pyx — coust . (47) 

From the definition of the pressure tensor we obtain 

r = _!ZM0i/2^, (48) 

2 dy 

p = ^i^[l + {l + e)Miy)]9. (49) 

Thus, we obtain the expressions for 9 and dux/dy as functions of p, r and v. Substituting them into the last equation 
of 14()|) we obtain 

A[^(.)|].G(.), (50) 

where 

n-) = v472[(J+-#(-'5))44-A(.)], (51) 

<^^) = e n 1-e , ,.,„ 52) 

with e = (t/p)2 and a{v) = 2/(i^[l + (1 + e)vg{v)]). 

It is well established to solve the second order ordinary differential equation such as H5U|) . Introducing H{u) as 
dH{v)/dv = F{v) and the multiplying dH/dy in both sides of (|5(J|>. and thus integrate the equation from y = to y 
we obtain 

where we use the symmetric condition dv / dy = dH/dy = at y = 0. 

± / ^ ' du = y. (54) 

^^0) ^2j;^^^^Fiv')Giu')du' 

Thus, we obtain the equation of y as the function of u. 

To draw the actual profile of v, we start from a trial J^i(O) to integrate (|54|l and calculate /i — /_^ /j ^i{y)dy, where 
the suffix 1 represents the first trial function. Then we replace vi{Q) by i^2(0) to reduce the deviation between /i and 
u. We repeat this relaxation procedure to obtain the converged result Im -^ v until M th trial. Once we obtain ;/, 
we can determine 9 and du^/dy by ea. (|48|) . 



13 




-90 -60 -30 30 60 90 

y 

(a) 





FIG. 9: The comparison between theory (sohd hnes) and DEM simulation without the tangential interaction (open circles) for 
area fraction (a), granular temperature (b) and Ux (c). The mean area fraction is i/ = 0.121. 
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FIG. 10: The comparison between theory (solid lines) and DEM simulation with the tangential interaction (open circles) for 
area fraction (a), granular temperature (b) and Ux (c). Here we use e — 0.798 and the mean area fraction !> = 0.121. 

To compare v and dux/dy with the result of DEM we use a fitting parameter e = (t/p)"^ , while we need two fitting 
parameters (t = —0.0017 and p = 0.1 for non-rotational cases and r = —0.0017 and p = 0.06 for rotational cases) to 
determine u^ and 9. It should be noted that p and t are determined by the boundary condition, but the boundary 
condition in eas. H29() and H31() with (|32|l contains two undetermined parameters. Thus, p and r cannot be determined 
within our theory. Figure |51 is the comparison of our theoretical result with the result of DEM without the tangential 
contact force in the interaction between particles for D = 0.121. The agreement between the theory and DEM is good. 
FigurelTUlis the comparison of our theoretical result with the result of DEM including the tangential interaction taking 
into account the renormalization of the restitution constant for D = 0.121. We again obtain a fairly good agreement 
between the theory and the simulation. 

The reason why we can use the kinetic theory is that collisions between particles are almost binary even in the 
dense cluster. Actually, we find that contacted particles are about 1.014 % of all particles at an instance, and only 
2.4 % of contacts are multi-body contacts among all contacts for D = 0.121. Therefore, the kinetic theory can be used 
even in the dense cluster in which particles are almost motionless. 

In principle, we can measure both the normal stress and the shear stress from the data of DEM. However, we only 
obtain the results with large errors. It seems that there is a tendency to have too small p in the direct measurement, 
though observed r is similar to the fitting value. 

Thus, we confirm that (i) the kinetic theory by Jenkins and Richman|l2|| can reproduce the profiles of hydrodynamic 
variables to describe the steady state of the granular fluid though the fitting values of the stresses are included, and 
(ii) the renormalization scheme proposed by Yoon and Jenkins[l3| is accurate. Although the setup of our simulation 
seem to be similar to that by Xu et al. 2!J| , our result for hydrodynamic variables is much heterogeneous than that by 
Xu et al. in the presence of a streamwisc flow. In fact, our system is separated into two regions which are a compact 
cluster and the very dilute region. The particles within the cluster is motionless, but the particles in the dilute region 
have large kinetic energy. On the other hand, the steady state obtained by Xu et al. is similar to USE. Thus, the 
result strongly depends on whether there is a stream flow. 



VI. DICUSSION 



Let us discuss our results in this section. To clarify the points we discuss three important points; the linear stability 
analysis of USE, the validity of the truncation at Navier-Stokes order and the effects of external fields. 
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A. The linear stability analysis 

Our analysis presented here suggests that USF of granular fluid is unstable. To verify such suggestion, we need to 
discuss the linear stability of USF. 

The stability of sheared granular flows has been discussed by many researchers. [Sll [S^, III, 154. ISSl 15 9 . 1571 IsM Is^ 
Many of them assume that the granular fluid is confined in an infinitely large system, but Alam et al. [S^] discuss the 
bifurcation of the steady solution as a function of the system size. We should note that the stability of granular flows 
on an inchned slope has also been discussed by some researchers. |6(1 iBlj 

Here, we have also checked the linear stability of the uniform sheared state assuming Lees-Edwards sheared boundary 
condition. Since the result is almost equivalent to that by Alam and Nott 53] , we only summarize the outline of our 
linear stability analysis. Our result indicates that USF is always unstable for enough large systems, but can be stable 
if the system size is enough small near e = 1. This result is consistent with the results by Babic[2^ and Popken and 
Clearv|24|. The details of our linear stability analysis will be reported elsewhere. 

As indicated in Section I, Bagnold's scaling can be used if the system is uniform and the heat conduction is negligible. 
However, the granular gas under the plane shear is not the case that we can assume Bagnold's scaling. The heat 
conduction plays an important role and USF cannot be maintained even when we adopt Lees-Edwards boundary 
condition. 

B. The validity of the approximation at Navier-Stokes order 

In this paper, we use the set of hydrodynamic equations at Navier-Stokes order. It should be noted that the set 
of hydrodynamic equations at Navier-Stokes order does not mean Newtonian. Actually, our granular system exhibits 
non-Newtonian behavior i.e. non-uniform velocity gradient as in Figs. 9 (c) and 10(c) in spite of the uniform shear 
stress in the steady state. 

From our analysis, the effects from the contraction of phase volume in collisions in the inelastic Boltzmann-Enskog 
equation arc also small. The effects of tangential contact force and the rotation of particles are also not important 
in the bulk behavior of hydrodynamics. Therefore, the kinetic theory by Jenkins and Richman 12] gives us sufficient 
accurate results to describe the hydrodynamics. 

Santos et al.^ and Tij et aL[62| indicate that the transport coefficients in Couette flow depend on the dimensionless 
shear rate. However, in our system the shear rate only determines the time scale and thus the qualitative behavior 
should not depend on the shear rate. There is room for discussion on the role of the shear rate in granular gases as 
an open problem. 

Some persons suggest the relevant role of terms at Burnett or super-Burnett order terms [SSL \63L l64| . but our analysis 
indicate that the contribution from higher oder terms should be small. It is known that hydrodynamic equations at 
Burnett order derived from the classical gas kinetics is unstable for perturbation, i.e. some solutions will blow up 
when a perturbation is applied to the system. [63, |6g| Therefore, we may need to be careful to use hydrodynamic 
equations at Burnett order or super-Burnett order to describe granular fluids. 

C. The effects of external fields 

Our system does not contain the external force except for the driving force of the boundary plates. In physical 
situations, it is difficult to remove the effect of gravity and collisions among particles may not be regarded as binary. 
As indicated in Introduction, recent papers and cited thereinjlfl llj of sheared granular materials under a constant 
pressure produce a state of 'granular liquid' in which particles are multiple contacts with each other. In these cases, 
the square root of stiffness of grains becomes a characteristic frequency. Thus, the behavior should depend on the 
shear rate. 

If we are interested in jamming transition or related phenomena of granular materials, we need to apply an external 
field to the system. Such subject will be important even in practical sense. 

VII. CONCLUSION 

In this paper, we have confirmed the validity of hydrodynamic equations derived from the kinetic theory by Jenkins 
and Richman. 12] We also confirm the relevancy of the renormalization method of the restitution constant by Yoon 
and Jenkins. [l3J This result may be surprised because the system includes a dense cluster whose packing fraction is 
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close to the maximum value, and the particles in the cluster are almost motionless. Since USF is unstable, we cannot 
use Bagnold's scaling to characterize the granular fluid in this case. 

We thank N. Mitarai and T. Hatano for fruitful discussion and their variable comments. HH also thanks J. T. 
Jenkins for his useful comments. This study is partially supported by the Grants-in-Aid of Japan Space Forum, 
and Ministry of Education, Culture, Sports, Science and Technology (MEXT), Japan (Grant No. 15540393 and No. 
18540371) and the Grant-in- Aid for the 21st century COE "Center for Diversity and Universality in Physics" from 
MEXT, Japan. 

APPENDIX A: THE DERIVATION OF WALTON'S /?o 

In this Appendix, we demonstrate how to derive Pq in Walton's expression in eq.|7I) 15] for the tangential restitution 
coefficient. Although the theory by Maw et a/.;16., 18, 67 , 6SJ has been used to evaluate /3, their expression is 
complicated and has an implicit form. Therefore, an explicit expression for Walton's /3o and the critical angle 7c is 
useful. 

Let us consider a collision between identical disks. Following the notation in Section II (without suffices i and j for 
colliding particles), the equation of motion for the tangential direction is described by 

wt + 2(r]twt + ktwt) = 0, (Al) 

when there is no slip during the collision. The factor 2 appears because the reduced mass is the half of mass of each 
particle. The solution of ea. (jAl|l is easily obtained as 

wt = ^^e-"'* sin(6i), mit) = u't(0)e-'"*(cos(6t) - ^ sin(6i)) (A2) 



for rif < 2kt, where b = \j2,kt — rj^- 

Since we choose large kt and small rjt-, the relation ry^ < 2kt should be satisfied. Similarly, the equation of motion in 
the normal direction is also described by an equation for a dumped oscillation. Thus, the duration time td at which 
Wn = is satisfied and the normal restitution constant are respectively given by 

^d^ /o. T' e = exp[ ]■ (A3) 

On the other hand, the restitution constant /3o for the tangential contact is thus given by 

Wtitd) 



/3o = 



miO) 



TrAT^t ,, 1 . ,TrATHy/2QjA^ 

exp( , ) , sm( , — ) 

77„V2i?^T'V2QM-l VnV2B~T 



where 



nArj^V^Q/A- 


')] 


' nnV2R - 1 


;Ji 







(A4) 



Trfe 



If we substitute the values fc„ — 3.0 x 10'^, kt = fcn/4, rjn — 3.0 and rjt — ?7„/2 used in DEM simulation, we obtain 
/3o ~ 0.769235. The comparison between the theory IJZJ) with ljA4|l and DEM is shown in Fig. El Without the 
introduction of any fitting parameters, agreement between the theory and DEM is fairly good. Here the critical angle 
cot 7c = (l + /3o)/i(l + e) ~ 1.56734. 
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